Dust Distribution in Gas Disks II: 
Self Induced Ring Formation Through a Clumping Instability 



Hubert Klahr 1 and D.N.C. Lin 
\J CO /Lick Observatory, University of California, Santa Cruz, CA 95064 
VERSION: Wednesday, June. 22nd 2005 
STATUS : ApJ in press 

ABSTRACT 

Debris rings of dust are found around young luminous stars such as HR4796A 
and HD141569. Some of these entities have sharp edges and gaps which have been 
interpreted as evidence for the presence of shepherding and embedded planets. 
Here we show that gaps and sharp edges in the debris disks of dust can also be 
spontaneously self generated if they are embedded in optically thin regions of 
gaseous disks. This clumping instability arises in regions where an enhancement 
in the dust density leads to local gas temperature and pressure increases. Con- 
sequently, the relative motion between the gas and the dust is modified. The 
subsequent hydrodynamic drag on the dust particles leads to further enhance- 
ment of their concentration. We show that this process is linearly unstable and 
leads to the formation of ring-like structures within the estimated life time of such 
young objects. Once the gas is removed (e.g. by photo evaporation) the struc- 
tures are "frozen" and will persist, even when the gas might not be observable 
anymore. 

Subject headings: circumstellar matter - planetary systems - stars: formation 



1. Introduction 

The coplanar orbits of nearly all planets in the Solar System inspired Laplace to propose 
that they are formed in a gaseous disk, the primordial solar nebula, centered around the 
protosun. Protoplanetary disks are found around a large fraction of classical and weak-line 
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T Tauri stars (Stauffer, Prosser, Hartmann, & McCaughrean 1994). The existence of these 
disks is generally inferred from the excess in the infrared continuum which is primarily due 
to dust reprocessing of the stellar radiation. These signatures fade with the age of the stars. 
In young stellar clusters, the fraction of stars with clear signs of infrared excess declines with 
their age on the time scale of a few Myr (Haisch, Lada, & Lada 2001). Around some post 
T Tauri stars such as (5 Pic, debris disks have been found. The amount of dust which emits 
infrared radiation in these debris disks also appears to be inversely proportional to the age 
of their host stars (Zuckerman, Kim, & Liu 1995). This apparent evolutionary pattern may 
be interpreted as coagulation of dusty material into relatively large particles (McCabe et al. 
2003) which are less efficient in reprocessing the stellar radiation (Yorke & Sonnhalter 2002). 
It is also possible that the residual grains are depleted from the circumstellar disks by photo 
evaporation, stellar wind, or some dynamical processes. In either case, the decrease in the 
infrared excess casts a constraint on the magnitude of time scale for planet formation to be 
in the range of a few million years. 

Current searches for extra solar planets indicate that they are present in at least eight 
percent of all target stars (Marcy, Butler, & Vogt 2000). As data accumulate with increasing 
time bases, the discoveries of long period planets are anticipated around a much larger 
fraction of all planets (Trilling, Lunine, & Benz 2002, Armitage et al. 2002, Ida & Lin 2003). 
If so, some debris disks may indeed contain newly formed protoplanets since they represent 
advanced evolutionary stages of protostellar disks. These planets are expected to tidally 
interact with their ambient gas and dust by exciting density waves (Goldreich & Tremaine 
1982, Papaloizou & Lin 1984), inducing gaps with sharp edges (Borderies, Goldreich & 
Tremaine 1989, Lin & Papaloizou 1986, Artymowicz, & Lubow 1994, Bryden, Ro'zyczka, 
Lin, & Bodenheimer 2000), confining rings between them (Bryden, Chen, Lin, Nelson, & 
Papaloizou 1999), and warping the disk (Larwood, & Papaloizou 1997). 

The discovery of a debris disk around HR4796A was particularly interesting because 
its sharp edges resemble that expected for tidally truncated disks (Schneider et al. 1999). 
Recent HST images of another disk around HD141569 clearly indicate that its outer region 
is strongly perturbed by its binary companion (Clampin, Ford, Illingworth, Krist & Ardila 
in preparation). But, this disk appears to be composed of two concentric rings which are 
separated by a narrow gap (Weinberger et al. 1999). The binary star's tidal torques occur 
on length scales comparable to the disk radius and are unlikely to induce the formation 
of the fine structure seen in this disk. A more plausible and potential exciting alternative 
culprit is a Jupiter-mass planet which may be embedded within the gap (Weinberger et al. 
1999). But, the distance of the hypothetical planets from their host stars in both cases are 
more than an order of magnitude larger than that of the gaseous giant planets in the solar 
system. The dynamical time scale in these disk regions exceeds 10 3 yr. Also in these regions, 
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the disk may attain sufficient mass to become gravitationally unstable and remain optically 
thin. The formation of giant planets through gravitational instability may be rapid (Boss 
1998) if the fragment can contract and collapse before being sheared out (Pickett, Mejia, 
Durisen, Cassen, Berry, & Link in press; Pickett, Durisen, Cassen, & Mejia 2000; Nelson, 
Benz, Adams, & Arnett 1998). The necessary condition for a rapid contraction is that the 
cooling time scale must be comparable to or shorter than the local dynamical time scale 
(Johnson & Gammie 2003). 

Application of the conventional core-accretion theory of planet formation suggests that 
the formation time scale for Jovian-mass planets at these large distances from the central 
object is much longer than the age of HR4796A and HD 141569 even under the most fa- 
vorable conditions (Pollack et al 1996). But the discovery of short-period planets (Mayor 
& Queloz 1995) supports the scenario that planets may migrate over extensive distance 
shortly after their formation (Goldreich & Tremaine 1980, Ward 1986, Lin & Papaloizou 
1986, Lin, Bodenheimer, & Richardson 1996). A particularly effective mechanism for rapid 
and extensive migration is through a runaway interaction between embedded planets and 
their nascent (Masset & Papaloizou 2003, Artymowicz in preparation) although the termina- 
tion conditions remain uncertain. Another potential avenue for extensive outward migration 
for protoplanets formed within a few AU from their host stars is through distant resonant 
encounter between several massive planets (Lin & Ida 1997). In the context of the Solar 
system, Thommes, Duncan, & Levison (2002) suggested that Uranus and Neptune may 
have migrated outwards as a consequence of the gravitational interaction between the outer 
planets and residual planetesimals. In general, dynamical instability induces the planets to 
attain large eccentricities. The gap width induced on the nearby disk by a highly eccentric 
planetary perturber may be much larger (Artymowicz & Lubow 1996) than that observed 
between the rings around HD 141569. But, the eccentricity of the perturber may also be 
rapidly damped by the collective response of the residual planetesimal disk (Agnor & Ward 
2002). 

The above discussions indicate that despite its potentially profound implication on the 
theory of planet formation, the planet-disk interaction scenario is highly uncertain. Thus, 
it is useful to explore some alternative explanations for the observed features in the debris 
disks. In a previous paper (Klahr & Lin 2001), we suggested that the debris disks and 
rings may be confined with relatively sharp edges through their interaction with both the 
stellar radiation and residual gas (sell also Takeuchi & Artymowicz 2001). Due to a negative 
pressure gradient in the radial direction, the azimuthal velocity of gaseous disks is generally 
sub Keplerian everywhere. In opaque regions of the disk, the grains' motion is affected only 
by the central star's gravity and their interaction with the ambient gas. As they attempt 
to attain Keplerian orbits, the grains experience a head wind and the hydrodynamic drag 
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resulting from this differential motion leads them to undergo orbital decay (Weidenschilling 
1977). 

But, in the optically-thin regions of the disk where the grains are directly exposed to the 
stellar radiation, the azimuthal velocity of the grains may also be reduced by the outwardly 
pointing radiation pressure from the central star. Around relatively luminous stars, such as 
HR4796A and HD 141569, the grains' azimuthal velocity may be less than that of the gas 
such that the combined effects of the radiation pressure and hydrodynamic drag can induce 
the grain to migrate outward. The above effects require the presence of residual gas which is 
otherwise difficult to detect directly. For example, at large distances from its central star, CO 
gas may be depleted as volatile molecules are condensed onto the grains (Thi et al. 2002). 
The detection of molecular hydrogen around some of these debris disks have been reported 
(Thi et al. 2001), although observational uncertainties remain (Sheret, Ramsey-Howat, & 
Dent 2002). 

In the context of the debris disk around HR 4796A, we showed that due to these hydro- 
dynamic drag and radiation pressure, a slight inhomogeneity in the gaseous background can 
lead to the formation of the dusty rings with sharp edges. Thus, the observed features in the 
ring do not necessarily imply the presence of embedded planets. Furthermore, if these sharp 
features are indeed an indicator for some residual gas, a less stringent formation time scale 
would be required for the formation of gaseous planets in the outskirts of planetary systems 
(for the formation of Neptune, see Bryden et al. 1999). 

Nevertheless, gas may have been present at the time when the structures in HR4796A 
and HD 141569 formed. After the ring structures were created by assistance of the gas, the 
gas may well have been evaporated and depleted to a level below observability and below a 
value that would still influence the dust. The dust rings would then be the "frozen" final 
state of the ring formation process. 

In this paper, we extend our previous investigation on the dynamical evolution of opti- 
cally thin debris disks around relatively luminous stars. In our previous work, we assumed 
an arbitrary surface density distribution for the gas. Our main objective here is to illustrate 
that the feedback interaction between the gas and dust naturally leads to a clumping insta- 
bility in which the dust drifts into an equilibrium configuration until the dusty rings become 
self confined. We show that this clumping instability results naturally and the structure of 
the debris disks can be determined self consistently. In §2, we outline our basic assumptions 
and formulate the evolution of both dust and gas. A linear stability analysis in §3 is followed 
by some numerical models in §4 and their observational signature in §5. We summarize our 
results and discuss their implications in §6. In the appendix we give a simple estimation for 
realistic f3 values (i.e. the sensitivity of gas temperature on grain density). 
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2. Basic assumptions for the model 



We consider an optically thin debris disk where the gas dust plasma is being irradiated 
by stellar light. The gas attains an equilibrium temperature somewhere between a maximal 
and a minimal possible temperature. The maximum temperature is given by the blackbody 
temperature of the dust grains and the minimum temperature by the ground state of the 
gas molecules and the temperature of the interstellar space. 



The equilibrium temperature of the gas T g in an optically thin regime at low temper- 
atures is given by the balance between radiative cooling via line emission and heating due 
to various effects including photo-electric heating via the silicate grains, the gas molecules' 
inelastic collisions with the dust grains as well as line absorption of UV photons (Kamp & 
van Zadelhoff 2001; Ferland, Korista, Verner, Ferguson, Kingdon & Verner 1998; Woitke, 
Kriiger & Sedelmayr 1996). The maximum temperature which can be attained by the gas 
is already given by the temperature of the dust grains, which themselves are heated to 
blackbody temperature by irradiation from the central star. In the current analysis, we are 
primarily interested in the interaction between the gas and the dust. We neglect effects such 
as super-heating of small grains due to radiative inefficiencies (Chiang & Goldreich 1997) for 
this analysis as it is not important to understand the basic instability. 

In general, the gas temperature is a complex function of radius R, luminosity of the 
central object L*, surface of the dust grains, surface density of the gas S 9 , surface density of 
the dust content E^, the cross section of the dust Ad, and most of all the complete chemical 
composition of the gas N n and the ionization state of the molecules X n , determining the 
possible line transitions leading to the emission of radiation: 



The detailed determination of the function / will be presented elsewhere. 

For the present purpose of examining the stability of dust to gas interaction, the most 
important parameter is 



In order to illustrate the evolution of small disturbances, we adopt a simple generic analytic 
prescription in which we assume that 



2.1. Gas Temperature 






(3) 
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where (3, T and S are to be determined later by radiation transfer models. In the linear 
stability analysis, only the instantaneous slope (3 needs to be specified and the above explicit 
prescription for T g (T, d ) is not needed. 

The power index (3 gives the steepness of the gas temperature dependence on the dust 
load. In regions where the dust to gas surface density ratio, Sd/Sc, > Z Q , where Z e is the 
solar metallicity, the disk gas is primarily heated by conduction to the gaseous molecules 
from the dust particles, which themselves are heated by the irradiation of their host stars 
and via photo-electric heating via the silicate grains. Qualitatively, T g is expected to be 
an increasing function of because in the optically thin region of the disk, the fraction of 
stellar radiation absorbed by the dust is proportional to the optical depth of the dust. The 
rate of conduction between the gas molecules and the dust grains as well as the number of 
photo-electric electrons is also an increasing function of S d . Thus, we consider primarily the 
range of parameter space for which (3 > (see appendix). However, in the low S^/E g limit, 
the gas molecules are mostly heated directly by the UV photons of their host stars such that 
the energy deposition rate to the disk is insensitive to the magnitude of (Kamp et al. 
2003). Also in the opaque regions of the disk, T g = Td and T g can be independent of 
so that (3 = 0. Therefore, we also consider other values of (3. Under unlikely circumstances 
where (3 < 0, the debris disks are stable against this clumping process. 



In our previous analysis (Klahr & Lin 2001), we derived the radial drift velocities v r for 
particles under the combined influence of non-Keplerian gas motion and radiation pressure. 
(The velocity of the dust is expressed in lower case v's whereas that of the gas is expressed in 
capital letters Vs.) For small particles with radius ad and density pd that are in the Epstein 
regime the friction time is 77 = ddPd/ (p g c s ) in a gas of density p g at a sound speed of c s . 
The radial velocity is then 



where the Stokes Number St = TfQ* « 1, Q* is the particles' hypothetical orbital frequency 
subjected only to the stars' gravity and radiation pressure. The difference of the actual 
orbital velocity of the disk gas and the hypothetical particle velocity Q*R, dV*, can have 
either positive or negative signs depending on the structure of the gas disk. 

Similarly, the radial drift velocity of the larger particles can be expressed as 



2.2. Radial Drift velocity 



v r = T f 2Q*dV* for St«l, 



(4) 



v r = dV* for St = 1, 



(5) 
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and 

v r = — for St » 1. (6) 

We utilize this expression for the radial drift velocities in the linear stability analysis below. 



3. Stability Analysis 

The evolution of an axis-symmetric dust distribution is determined by the continuity 
equation in cylindrical coordinates: 

d 1 d 

-£, + -— r£^ r = 0. (7) 
at r or 

In general one would have to solve a continuity equation for all single dust species as their 
radial drift velocity depends on the friction time. Such detailed approach will be carried out 
in future models which will try to fit the observations in detail. For this work we assume 
a mono dispersive dust distribution which is characterized by one particle size, shape and 
friction time. 

The radial drift velocity for the most mobile particles (St = 1) is: 

v r = dV\ (8) 

and dV* — — is a function of the particles theoretical azimuthal velocity for circular 
orbit in the absence of gas 

4nm d c /g\ 

r 

Here A d is the cross section and rrid the mass of the dust grains. The speed of light is denoted 
by c. 

We want to stress that the characteristic radial drift velocity will also be correct within 
an order of magnitude for particles which are up to 20 times larger or 20 times smaller than 
the characteristic size. 

The azimuthal velocity of the gas is and is determined by the radial pressure gradient 
in first oder approximation to 

where f2 = (GM./r 3 ) 1 ' 2 is the Keplerian frequency. The pressure P is given by an ideal gas 
equation of state a function of temperature and gas surface density S s thus: 
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where R gas is the gas constant. Note that the second and third terms on the right hand side of 
eq. (11) are due to the surface density and temperature distribution of the gas respectively. In 
general, £ 9 is not perturbed by the dust concentration. But any clumping in the dust would 
lead to a local enhancement in which in term heats the gas. Because the temperature 
distribution is much more sensitive than the S s distribution, we focus our stability analysis 
on those contributions which are proportional to the temperature gradient. 

Combining all the equations into one expression and defining a local mean sound speed 
to be ^ill = c 2 we fi n( l; 



<9 V _ Id 

777 — 77-?" 

ot r or 



2f2r J \T J V dlnr dlnr 



• (12) 



With a prescription for T g such as that in eq(3), this continuity equation can be solved 
numerically. 

For the linear stability analysis, we simplify the equation to clarify the mechanism by 
dropping all terms not important for the instability. We neglect the effects of radiation 
pressure which gives only a systematic offset but no instability. We also assume that the gas 
is much more smoothly distributed than the dust, i.e. |<91n£ 9 /<91nr| << \d\nT, d /dr\ so that: 

d ^ 1 d f (T g \ c 2 s (3d\nE d 



dt^ rdr\\T ) 2VL d\nr )' ^ 

Linearized = S + S^, we find in a local approximation with ^d r r — > dx at a given 
radius r = R 

d c 2 d 2 

m^ + ^^ = °- (14 > 

The local approximation does not limit the physical validity of our analysis. It simplifies the 
mathematics but does not change the underlying physics. Using a polar coordinate system 
would have required to apply e.g. Bessel functions instead of Fourier modes. The result 
would have been the same but it would be less transparent. However, {l/rd r r = d x ) requires 
x/r < 1. In the case of HR4796A the radial extent of the ring is 20 AU where the radius is 
70 AU, so the approximation is fulfilled in this system. 

For short- wavelength disturbances, the above equation has a solution S' d ~ e -iM-fcx)_ 
The growth-rate of these disturbance is: 

r = -iu = (3^k 2 . (15) 
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If we replace c s = HQ, and define n = kR then the above equation reduces to 

r = - <w = ^( n f) a ( 16 ) 

This result indicates that a S 9 variation can lead to local heating of the gas which in term 
modifies the local pressure gradient and the azimuthal speed of the gas. Through a hydrody- 
namic drag, the dust particles respond by clumping together to enhance the local S fl . This 
clumping instability is related by, but not the same as the diffusion and thermal instability 
in gaseous accretion disks (Pringle 1981). In both cases, the growth of a perturbation is due 
to the efficiency of angular momentum transfer being inversely proportional to the surface 
density. But the present clumping instability involves a different physical process which re- 
quires the feedback between the gas and the dust through the radiative heating of the dust 
and the local heating of the gas via the dust. 

It is interesting to notice that the growth rates are close to the orbital frequency for 
H ~ K In this limit, the instability will grow on a dynamical time scale. For disturbance 
with n > R/H the linear leading-terms approximation is inappropriate. On these scale, the 
thin disk approximation breaks down and the gas pressure in the direction normal to the 
plane has a tendency to stabilize the flow. In addition, very small scale variations, if unstable 
can reverse the gradient of angular momentum in the disk (see eq. 11). A negative angular 
momentum gradient, even localized can lead to interchange instabilities which results in 
mixing in the radial direction. Since we are primarily interested in identifying the condition 
for the onset of a clumping instability, we shall present elsewhere the growth limit due to 
amplitude saturation and nonlinear effects. 

Nevertheless, eventually the gas will be removed from the system (e.g. photo evapora- 
tion). As a consequence the friction time of the particles will increase and the radial drift 
velocity decrease (see eq. 6). Thus, the right hand side of eq(13) will eventually vanish. This 
means that the structure is now in first order approximation constant in time. So even the 
gas is removed from the system the structure in the dust distribution will remain. Conse- 
quently we can not argue that there must still be gas at a dynamically important or even at 
an observable amount in HR4796A and HD141569. 

Note that T and (3 have the same sign. It is not forbidden that (3 may attain a negative 
value, in which case the perturbation would be damped. However, this case requires the gas 
temperature to decrease with which seems to be physically unlikely for an optically thin 
debris disk. Only when the dust grains are the efficient cooling agent and the heating of the 
gas is provided by other means then f3 might become negative. Yet this case is unlikely for 
the objects that we consider. 

In opaque regions of a circumstellar disk where all radiation is intercepted by the dust, 
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the gas temperature becomes independent of so that (3 — 0. In these (3 — cases, the 
perturbation remains neutral and the disk is stable. 



4. Numerical Solution 

In order to have an illustration of how the instability can evolve into the non linear 
region, we obtain numerical solutions of the complete continuity equation 12 without sub- 
stituting the gas temperature T g and assuming circular orbits: 

Here we set all physical constants to one and used VL = r~ 3 / 2 . Using these dimensionless 
units sets the unit time scale to one orbital period at the radius r times the inverse square 
of the pressure scale height. 

t "=\t) a < 18 > 

Thus, it will be simple to calculate the physical time scale for a real object like HR4796A 
afterwards by simply multiplying the dimensionless time by the unit time scale of the physical 
radius, that is considered. 

The gas temperature is calculated from the power-law prescription for T g as expressed 
in eq. (3), where we additionally limited the dynamical range to one order of magnitude: 

3.33, min 



T g = T max 



0.33, ( g 



(19) 



As we will discuss in the appendix the gas temperature will probably stay in such limits (see 
Figure 1). 

The analytic estimate in eq. (16) indicates that the growth rate is proportional to n 2 , 
i.e. the most unstable perturbations are those with the smallest wave numbers. We have 
already indicated above that the thin disk approximation for the continuity equation breaks 
down for n > R/H. In reality, the vertical pressure gradient becomes important to stabilize 
the disturbance. We also argued that the non linear growth of disturbance with n > R/H 
would lead to the violation of the Rayleigh's criterion for axis-symmetric rotational flows. 
If such a situation arises, mixing would almost certainly limit the growth of the clumpy 
structure. In order to avoid the unphysical growth on scales smaller than H, we introduce a 
prescription by adding a stabilizing contribution to the continuity equation such that 
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With this prescription, the equation can now be discretized with central differences 2 on a 
grid with 100 radial points and 2 ghost cells on each side, where T, d was set to fixed values 

dinner = 1-0 and Yi ou i er = 10 10 . 

Another choice of boundary conditions would have been vanishing gradients. But this 
leads to numerical problems with the scheme we applied, i.e. negative densities can occur 
close to the boundary. The application of fixed values for the ghost cells helps to avoid this, 
whereas the overall ring structures are not changed in size and location. 

Once one starts to model a given observation, one has to include more physics than 
we do now. Then one needs the detailed cooling and heating behavior for the optical thin 
gas. One will then also find that only a certain region around the star will be subject to 
the clumping instability. For instance at very large radii the gas will be at a lowest possible 
temperature of maybe 10K which will be determined by the back ground radiation field, thus 
there is no feed back mechanism from the dust density to the temperature of the gas. At radii 
to close to the star dust might be destroyed by high temperatures and also no such instability 
is possible. In reality the radial window for instability might be even narrower than that, 
because the cooling behavior depends on the chemical composition of the disk gas, which 
itself depends on the radiation field. Only where the right cooling agents are present one can 
hope that the cooling behavior gives rise to the clumping instability. Detailed investigation 
on real objects are necessary to identify the potentially ring forming regions. 

In such a system, where the instability is truncated towards the edges, one will have 
more natural boundary conditions. In our simplified model we thus assume that the dust 
density does not change at the edges, which has a similar effect as shutting of the instability 
there. 

We use a numerical scheme using a fourth order Runge Kutta integrator to advance 
the distribution in time. The fixed background gas density distribution was chosen to be 
Tjg oc r~ 3//2 . A similar profile was chosen for initial distribution of the dust with a 
randomly fluctuation of 1% amplitude. We also set T oc r -1 / 2 with H = 0.1. 

The stabilizing contribution is the only limit for the growth of the shortest modes in 
this flat approximation where vertical structure is neglected. A more realistic modeling 
would require to resolve the 3D structure of the gas disk and the turbulent processes within. 
Nevertheless, this is way beyond the considerations of this paper. On top of this we would 



2 The necessity for the stabilizing contribution arises from the fact that we use a central differencing 
scheme. In future projects we will adapt an upwind advection scheme which is stable by construction even 
without the artificial stabilization. 
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only shift the application of free parameters one layer up, without gaining much additional 
insight. 

A pressure change can lead to the modification in the density distribution in both vertical 
and radial direction. Since hydrostatic equilibrium is maintained in the vertical direction 
by a balance between pressure and gravity, a change in the pressure can modify the density 
distribution above the mid plane. But it does not necessarily change the surface density. 
In the radial direction, equilibrium is maintained by a balance between centrifugal force of 
the Keplerian rotation and the stellar gravity. Thus, small variations in the pressure do not 
modify the gas density significantly. 

The linearized solution in eq(16) clearly indicates that the sign of /3 determines the 
growth and damping of the disturbances. We choose three values of (3 (-1, 1, 2) to illus- 
trate this dependence. The magnitude of f3 also modifies the radial location of the dust 
concentration maxima. 

Fig. 2 shows the evolution of a small random fluctuation with (3 = 1. The solid line 
denotes the final time, where we stopped the integration at the dimensionless time t max = 1. 
The density is only fixed in the ghost zones, which we do include in the plot, thus the values 
close to the boundary may differ from the fixed values (see above). The doted lines are 
snapshots after evenly spaced periods of time. We plotted the surface density of the dust 
S rf (t) over the dimensionless radius. 

The instability is such an effective process that the initial state does not matter. With 
more effort we could have started with a smoother distribution towards the edges, but the 
result would have been the same. 

The results in fig. 2 clearly indicate the growth of the perturbed and the formation of 
ringlets which are separated by gaps. The width of the rings is an increasing function of the 
radius because H/r increases with r. Eventually, as the amplitudes of S d are strong enough, 
the residual disk attains an equilibrium in which the residual dust rings become isolated. 

The numerical models needed 1 dimensionless time scale to develop rings. The mean 
radius of the simulation was chosen to be r = 5 thus it took only 0.1 dimensionless orbits at 
the center of the dust density ring. This typical time scale can be derived from eqs. (12) and 
(19) to be: r = (— ) 2 ^y. If we apply parameters estimated for HR4796A (see Model A in 
Klahr and Lin 2001) we can estimate the time the instability needed for growth at the ring 
location of 70AU. The orbital period is 370 yrs. The pressure scale height is determined from 
the irradiation to be H/R = 0.07. It follows that the typical time scale is only 7.5 x 10 4 yrs, 
which means the instability developed in the model after about 10 4 yrs for the most efficient 
drifting particles. This model (Model A from Klahr and Lin 2003) assumes the existence of 
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gas at an amount of 100M e in the disk around HR4796A. With such a high gas density dust 
grains of 600yU in size are drifting radially at the highest possible rate (fir/ ~ 1 see eq. 8). 
The radial drift velocity decreases with the friction time and thus with gas density. In Klahr 
and Lin (2001) we also assumed two models (Models C and D) that contain only 10 M® 
respectively 1 M® in the gas. Under such conditions the instability needs 10 respectively 100 
times longer. For HR4796A this means that the ring may have formed within 10 4 and 10 6 
yrs depending on the amount of gas that was available. Thus, if we assume that the gas was 
depleted on a time scale of 10 6 yrs it is possible to generate the ring structure within this 
time and also to conserve or fixate the structure by removing the gas during the process. 

In Figure 3 we present the evolution of for the (3 — 2 model. Maximum time was also 
tmax = 1- The growth of the perturbation in this case is similar to the (3 = 1 model though 
the growth rate is faster and the outward drift of the rings is more pronounced. Finally we 
present in Figure 4, the evolution of S d for the (3 = —1 model. As this process is generally 
slower we integrated until t max = 10. As we have anticipated in our analytic results that an 
initially strong sinusoidal perturbation is damped in this case. 

5. Observational Appearance 

The radial dust density and temperature distribution can be translated in an intensity 
map via 

7(r,0)ocE d T d 4 . (21) 

In figures 5 and 6 we show these maps for the (3 = 1 and (3 = 2 case. In the latter case 
the middle ring is a little wider than in the previous case. One should not expect to much 
difference between both cases, as (3 predominantly influences the time scale on which the 
rings form but not the structures themselves. 

The resemblance to the observations of HR4796A and HD141569 is striking. Future 
work could focus on a possible instability of non-axisymmetric modes, which are actually 
already observed. 

6. Conclusion 

In this paper, we consider the interaction between dust, gas, and radiation in debris 
disks around young stellar objects. In general, the gas and the dust component do not 
have the same azimuthal speed. The gas attains a dynamical equilibrium in which the 
centrifugal force is balanced by the central stars' gravity and a negative pressure gradient in 
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the disk. The azimuthal speed of the gas is generally sub Keplerian. In a gas-free optically- 
thin environment, the dust attains a dynamical equilibrium in which the centrifugal force is 
balanced by the central stars' gravity and a radiation pressure. In a gaseous environment, 
hydrodynamical drag also influences the motion of the gas and dust. In most regions of 
the disk where the dust concentration is relatively small, the gas causes the dust to migrate 
inward/outward if the azimuthal speed of the gas is smaller /larger than that of the dust. 

However, in relatively massive Herbig Ae/Be stars, the UV intensity is relatively strong. 
Gas may be rapidly depleted through photo evaporation process (Shu, Johnstone, & Hollen- 
bach 1993, Hollenbach et al. 1994). Nevertheless, a small amount of residual gas and dust 
may be retained around systems such as HR4796A and HD 141569 for a few Myr. Regions 
with modest gas density and dust opacity are primarily heated by the absorption of stellar 
irradiation on the dust grains. The gas is indirectly heated through photo-electric heating 
via the dust grains. 

We consider the evolutionary stability of a two-fluid flow under these circumstances. 
For simplicity, we consider the limiting case that the azimuthal speed of the gas is smaller 
than that of the dust initially. (Similar results can be obtained in the opposite limit). We 
find that in some optically thin regions of the disk, the gas temperature is an increasing 
function of the surface density of the dust. In this case, a slight concentration of dust leads 
to local pressure maxima in the gas. The positive pressure gradient interior to these maxima 
increases the local azimuthal velocity of the gas which remains smaller than that of the 
dust initially. Although the dust continues to experience a hydrodynamic drag and migrate 
inward, the strength of the head-wind is weaker and the inward radial velocity is smaller. 
Similarly, exterior to these maxima, the pressure gradient is more negative than elsewhere 
so that the azimuthal speed of the gas is slower than its surrounding region. This enhanced 
head wind causes the dust to migrate toward the maxima at a faster pace. The congestion 
caused by the slow down in the radial migration of the dust on the inside and the speed up 
on the outside of the surface density maxima in turn leads to an increase in S^. This process 
leads to a clumping instability. 

The initial growth of the clumping instability is similar to the viscous diffusion instability 
in accretion disk models. But, in contrast to the viscous instability, the rings can evolve into 
apparently isolated entities. As the amplitude of the surface density variation increases, the 
modification of the local pressure distribution intensifies. Eventually, the azimuthal speed 
of the gas decreases below that of the dust in some region interior to the surface density 
maximum. The dust particles move outward as they experience a tail wind. This inflow 
barrier leads to the apparent isolation of the rings, while the gas distribution is merely 
unchanged during the entire process, but the gaseous component of the residual disk is 
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difficult to be detected anyway. Although this self shepherding effect diminishes with gas 
depletion, the diffusion of the grains due to the gas drag effect also vanishes. After the 
residual gas is completely evaporated, the ring structure would be preserved until particle 
collisions and radiation drag induces its spreading. But these processes occurs on a much 
longer time scale than a few Myr. 

The growth rate of this clumping instability is inversely proportional to the radial wave- 
length of the perturbation so that it leads to the formation of concentric rings with gaps 
between them. The lower limit on the width of the ring is determined by the local pressure 
scale height. In accordance with the Rayleigh's criterion, rings with narrow widths are dy- 
namically unstable. Although our discussion is based on the assumption that the azimuthal 
velocity of the gas is initially smaller than that of the dust, the same argument applies for 
the opposite limit. In that case it is the congestion of the dust's outward flow which causes 
the growth of the clumps. 

Debris disks have been seen around many post T Tauri stars. Some of these systems 
such as (3 Pic have relatively smooth surface density distribution in the radial direction 
while those around HR4796A and HD141569 appears to have ringlet structure. While it is 
tempting to interpret the ringlet structure with a tidal interaction between the residual disk 
and some embedded planets, our results indicate the existence of alternative scenarios which 
can also account for the sharp features in the ring. 

The stability condition is determined by the functional dependence of T g on S^. The 
instability arises when the gas temperature is an increasing function of the dust surface 
density. Such conditions are not satisfied in the optically thick regions of the disk or in 
regions where the gas molecules radiate either much more or much less efficiently than they 
are heated by conduction. If gas is an efficient radiator, T g would attain the value that 
corresponds to the ground state of the gas molecules. In the opposite limit, gas would attain 
the radiative equilibrium temperature of the grains. In both limits, the gas temperature 
becomes independent of S d , (3 = 0, and the debris disk is stable. It is also possible that 
(3 may attain a negative value though it seems physically unlikely for a debris disk. We 
shall examine elsewhere the functional dependence of the gas temperature on other physical 
parameters. This analysis will include a more realistic dust size distribution and the radiation 
transfer of the gas. We will also provide detailed modeling of the observations. It is entirely 
possible that only over some special range of disk radii the clumping instability may arise. 

We want to thank Inga Kamp for her thoroughly checking of the heating and cooling 
properties of optically thin disks. This work was supported by NASA TPF program under 
the grant NNG04G191G. Work supported in part by the European Community's Human 
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Potential Programme under contract HPRN-CT-2002-00308, PLANETS. 

A. Gas Temperature Estimate 

For this work we assumed a temperature versus dust density relation with a positive 
slope, without further discussing the exact shape of this relation. A detailed exact derivation 
would need the use and modification of a radiation transfer code (see Kamp & van Zadel- 
hoff 2001; Ferland, Korista, Verner, Ferguson, Kingdon & Verner 1998; Woitke, Kruger & 
Sedelmayr 1996). This shall not be done here. 

But we can derive some general expressions for the gas temperature (T) as it depends 
on the surface area that is provided by the dust grains. The strongest heating Q + provided 
to the Hydrogen molecules probably occurs by photo-electric heating via the silicate grains 
(Kamp & van Zadelhoff 2001): 

r pe = 2.5 x 10 -4 <rexn H . (Al) 

Here a = is the total grain cross section per H nucleus, e is the photo-electric efficiency, 
and x is the intensity of the UV field from the central object (a^ is the cross section and n 
the number of dust grains). 

Actually e also depends on the gas temperature, but for the expected low temperatures 
T g « 10 4 K this dependence is negligible for our purpose, e depends on n e and x an d hence 
using reasonable values of those two, e.g. n e = 1cm -3 , x — 300 at 70 AU, e can only vary as 
much as by a factor 2 between 10 and 100 K. 

Thus it follows, that the heating is in first order linear increasing with the number of 
the dust grains: 

T pe oc n. (A2) 

The Hydrogen being at its ground state can not cool via radiation but only via collisions 
with other atoms or molecules. These collision partners themselves will cool via rotational 
or fine structure transitions (Woitke et al. 1996). Given the strong UV radiation field of a 
young A0V star, the dominant cooling agent is probably ionized carbon CII (Kamp et al. 
2003). One can expect that C is totally ionized and nc+ oc 10~ 4 nH- We use the two level 
approximation of Hollenbach & McKee (1989) for the ground state of the ionized carbon: It 
is the 157.5 \x line that has a cooling rate of: 



A 2 = /in(CII)Ai /w/io. 



(A3) 
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The fraction of CII in the upper level can be calculated via 

/! = — — , (A4) 

with the statistical weights E = 0.0 and E\ = 1.27 x 10~ 14 and go = 2.gi = 4 respectively. 
A w is the Einstein probability for spontaneous emission. When the heating is balanced by 
cooling, we can calculate the dependence of equilibrium temperature on the number of dust 
grains: 

9ie- E ^ /kT n 
go + gie-*/«r « - 

which can be easily solved for T. In Fig. 1 we plot the dependence of T/T Q on n/n . The 
exact values for T and n have to be calculated model dependent, but we already see that 
(3 in this plot is roughly 1.0 (which is independent on the explicit model) and we conclude 
that the instability can occur on reasonable time scales. 
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2.5 



Fig. 1. — Gas temperature in dependence on the number of dust grains. See text for details. 




Fig. 2. — Evolution of an initially homogeneous dust distribution over dimensionless time 
t — 1 in units of the initial value over the dimensionless radius. (3 — 1. 




Fig. 3. — Evolution of an initially homogeneous dust distribution over dimensionless time 
t — 1 in units of the initial value over the dimensionless radius. (5 — 2. 




Fig. 4. — Evolution of an initially sinusoidal dust distribution over dimensionless time t — 10 
in units of the initial value over the dimensionless radius. j3 — — 1. 
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Fig. 6. — Simulated intensity map for an optically thin debris disk. (3 — 2. Now the middle 
ring is wider than in the previous case. 



